1. Root Finding#
강좌: 수치해석 프로젝트
이론#
비선형 방정식의 해를 찾는 과정으로 반복적으로 해석한다.
대표적인 방법은 Bisection Method, Newton Rapson Method가 있다.
Bisection Method#
해가 구간 (a, b) 사이에 있을 때 해가 존재하면 \(f(a)\cdot f(b) < 0\) 조건을 활용한다.
이 경우 중간 값 \(c=\frac{a+b}{2}\) 에 대해
\(f(a)\cdot f(c) < 0\) 이면 해가 (a, c) 사이에 있다고 범위를 좁힘
\(f(c)\cdot f(b) < 0\) 이면 해가 (c, b) 사이에 있다고 범위를 좁힘
Note
\(|f(c)| < tol\) 이면 멈춘다.

Fig. 1 Bisection method (From Wikipedia)#
Python 구현#
def bisect(a, b, f, tol=1e-12):
"""
Bisection method
Parameters
----------
a : float
Lower limit
b : float
Upper limit
f : function
함수
tol : float
Tolerance
"""
product = f(a)*f(b)
if product > 0:
# 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
raise ValueError('Wrong Intervals')
else:
c = 0.5*(a+b)
if abs(f(c)) < tol:
print('Converged at {:.7g}'.format(c))
elif f(a)*f(c) < 0:
# 해는 (a,c) 사이에 존재함, 범위 좁히기
bisect(a, c, f, tol)
else:
# 해는 (c,a) 사이에 존재함, 범위 좁히기
bisect(c, b, f, tol)
예제#
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
x = np.arange(-10, 10, 0.1)
def f(x):
return x**2 + 10*np.sin(x)
plt.plot(x, f(x))
[<matplotlib.lines.Line2D at 0x7f49dc146210>]

# Left solution
bisect(-10, -1, f)
Converged at -2.479482
# Right solution
bisect(-1, 1, f)
Converged at 0
# Wrong Interval 1
bisect(-10, -5, f)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In [7], line 2
1 # Wrong Interval 1
----> 2 bisect(-10, -5, f)
Cell In [1], line 20, in bisect(a, b, f, tol)
16 product = f(a)*f(b)
18 if product > 0:
19 # 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
---> 20 raise ValueError('Wrong Intervals')
21 else:
22 c = 0.5*(a+b)
ValueError: Wrong Intervals
# Wrong Interval 2
bisect(-5, 5, f)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In [8], line 2
1 # Wrong Interval 2
----> 2 bisect(-5, 5, f)
Cell In [1], line 20, in bisect(a, b, f, tol)
16 product = f(a)*f(b)
18 if product > 0:
19 # 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
---> 20 raise ValueError('Wrong Intervals')
21 else:
22 c = 0.5*(a+b)
ValueError: Wrong Intervals
Newton Rapson Method#
함수 \(f(x)\) 가 미분 가능하고 연속함수인 경우에 대해서 다음과 같은 근사식을 만족한다.
즉,
이 과정을 반복하면 해를 구할 수 있다.
Note
\(|f(p_n)| < tol\) 이면 멈춘다.
초기 Guess 값에서 시작해야 하며, 이에 따라 해가 달라진다.
방정식의 미분 \(f'(x)\)가 필요하며, 이론적으로 구하기 어려운 경우 수치적으로 계산해야 한다.

Fig. 2 Newthon Method (from Wikipedia)#
Python 구현#
def newton(f, df, x0, tol=1e-9):
"""
Newton Rapson method
Parameters
----------
f : function
함수
df : function
도함수
x0 : float
초기 Guess
tol : float
Tolerance
"""
# Do It Yourself
if abs(f(x0)) < tol:
print('Converged at {:.7g}'.format(x0))
else:
# x 절편 (Intercept)를 계산함
x0 -= f(x0) / df(x0)
# 새로운 점에 대해 Newton Rapson 기법 적용
newton(f, df, x0, tol)
예제 해석#
도함수는 이론적으로 또는 수치적으로 계산할 수 있다.
def df_exact(x):
return 2*x + 10*np.cos(x)
def df_apprx(x, tol=1e-10):
return (f(x+tol) - f(x-tol)) / (2*tol)
# Left solution
newton(f, df_exact, x0=-4)
Converged at -2.479482
# Left solution with
newton(f, df_apprx, x0=-4)
Converged at -2.479482
# Right solution
newton(f, df_exact, x0=3)
Converged at 5.448369e-16
Scipy 활용#
scipy.optimize
모듈은 최소화, Curve fitting 그리고 root finding과 관련된 다양한 알고리듬을 제공한다.
참고
from scipy import optimize
scipy.optimize.root_scalar
또는 scipy.optimize.root
함수를 이용하면 손쉽게 계산할 수 있다.
Note
scipy.optimize.fsolve
함수는 MINPACK 라이브러리의 wrapper로 MATLAB의 fsolve와 용법이 같으나, 이제는 안 쓸 것을 권장한다.
# Root with Hybrid
optimize.root(f, x0=-4)
message: The solution converged.
success: True
status: 1
fun: [ 1.510e-14]
x: [-2.479e+00]
nfev: 8
fjac: [[-1.000e+00]]
r: [ 1.285e+01]
qtf: [ 1.069e-08]
# Root_scalar with brentq method
optimize.root_scalar(f, bracket=[-10, -1])
converged: True
flag: 'converged'
function_calls: 12
iterations: 11
root: -2.479481833541344
# Root_scalar with Newton Rapson method
optimize.root_scalar(f, x0=-4, fprime=df_apprx)
converged: True
flag: 'converged'
function_calls: 10
iterations: 5
root: -2.479481833541344
# Find minimum from x0=0
optimize.minimize(f, x0=0)
message: Optimization terminated successfully.
success: True
status: 0
fun: -7.945823375615215
x: [-1.306e+00]
nit: 5
jac: [-1.192e-06]
hess_inv: [[ 8.589e-02]]
nfev: 12
njev: 6
# Find minimum from x0=5
optimize.minimize(f, x0=5)
message: Optimization terminated successfully.
success: True
status: 0
fun: 8.31558557947746
x: [ 3.837e+00]
nit: 5
jac: [ 2.384e-07]
hess_inv: [[ 1.188e-01]]
nfev: 12
njev: 6